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ABSTRACT 

This  paper  is  concerned  with  the  computation  of  semiconductor  device 
current-voltage  characteristics.  We  describe  an  algorithm  which  allows  the 
computation  of  characteristics  by  continuation  in  a  parameter  which 
approximates  the  arc length  of  the  characteristic.  The  use  of  this 
parameterization  allows  the  characteristic  to  continue  beyond  snap-back- 
voltages,  while  continuation  in  the  voltage  fails  past  snap-back-voltages.  We 
discuss  the  implementation  of  the  parameterization  and  give  a  numerical 
example . 
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SIGNIFICANCE  AND  EXPLANATION 


\  s-\ 

Injthis  paper  we  present^an  algorithm  for  the  numerical  computation  of 
current-voltage-characteristics  in  a  semiconductor  device  (i.e.  curves 
describing  the  dependence  of  the  output  current  on  the  input  voltage  at  some 
contacts).  For  certain  devices  (e.g.  thyristors)  the  set  of  equations 
describing  the  state  of  the  semiconductor  device  have  multiple  solutions  for 
specific  voltage  ranges  and  snap-back  phenomena  occur  (i.e.  the  current  is  not 
a  single  valued  function  of  the  voltage).  This  implies  that  the 
characteristic  cannot  be  computed  numerically  by  merely  stepping  up  the 

TU 

voltage,  ©or  algorithm  involves  reparameterizing  the  characteristic  curve  by 
means  of  a  parameter  which  approximates  the  arclength  of  the  curve. 
Continuation  in  this  parameter  past  snap-back-points  does  not  cause  serious 
difficulty. 

_We-  discuss  the  implementation  of  the  reparameterization  and  present  a 


numerical  example. 


The  responsibility  for  the  wording  and  views  expressed  in  this  descriptive 
summary  lies  with  MRC,  and  not  with  the  authors  of  this  report. 


COMPUTATION  OP  CURRENT- VOLTAGE-CHARACTERISTICS  IN  A  SEMICONDUCTOR  DEVICE 
USING  ARC-LENGTH  CONTINUATION 

Peter  A.  Markowich* ' 1 ,  Christian  A.  Ringhofer1,  and  Alois  Steindl* 

1 .  Introduction 

The  coaputation  of  static  current-voltage  characteristics  for  devices  which  exhibit 

snap-back  phenoaena  (like  thyristors  or  even  a  pn- junction  in  the  avalanche  case)  has 

created  probleas  since  continuation  in  the  voltage  does  not  work  beyond  snap-voltages. 

Figure  1  shows  a  typical  characteristic  of  a  thyristor  in  forward  bias  with  the  snap-back 

voltages  U_  and  U-  . 

S1  s2 


Figure  1.  Current  Voltage  Characteristic  of  a  Thyristor. 
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Amum  that  tha  point  (U,,!,)  -as  computed  (say  by  starting  at  (0,0)  and  by 
gradually  stopping  up  tha  voltaga  by  tha  increment  AU>.  Than,  if  U2  -  U,  +  AU  >  U^, 
tha  point  (02.12)  alraady  lias  on  tha  uppar  branch  of  the  characteristic  and  Newton’s 
Mthod  for  tha  numerical  «*put.tion  of  (02,I2>  -ill  probably  fail  -h.n  providing  an 
initial  guass  which  is  computed  just  by  using  tha  previously  computed  solutions.  Tha 
(dynamically  stabl.)  upper  branch  — rging  from  cannot  be  computed  by 

continuation  in  O  starting  fro-  soma  point  on  the  (also  dynamically  stable)  ’lower 

branch*  connecting  (0,0)  and  (°Sj' IS2 ' ' 

In  this  paper  wa  demonstrate  ho-  to  avoid  this  problem  by  introducing  an  arclength- 
type  parameter.  Tha  characteristic  can  be  computed  past  snap-back  voltages  when  doing 
continuation  in  this  new  parameter.  This  procedure  is  well  known  to  mathematicians, 
however,  most  semiconductor-device  simulation  codes  do  not  make  use  of  it  yet. 

We  also  demonstrate  how  an  already  existing  code,  which  solve,  the  fundaments i 
semiconductor  device  equations  for  given  contact  voltages  has  to  be  augmented  to 
accommodate  the  reparam.fr  ir  at  ion.  The  augmentation  turns  out  to  be  very  cheap  in 

progressing  as  well  as  in  computer  resources. 

Finally  we  present  numerical  results  for  a  one-dimensional  pn-junction  in  the 

avalanche  case. 
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2.  State— nt  of  the  Problem 

The  basic  static  semiconductor  device  equations  as  given  by  Van  Roosbroeck  (19S0)  are 

(2.1)  div(eV»>  -  n  -  p  -  C 

(2.2)  div(D  Vn  -  k  n7f)  -R'  xeQCRn,  n  -  1,  2  or  3  . 

n  n 

(2.3)  div(DpVp  +  upV»)  -  R  . 

(He  employ  the  usual  notation)  where  0  represents  the  device  geometry.  The  electron  and 
hole  current  densities  are  given  by 

(2.4)  Jn  -  q(DnVn  -  MnnV*)  , 

(2.5)  Jp  -  -q(DpVp  +  UppV>>) 
and  the  total  current  density 

(2.6)  J  -  J_  +  J_  . 

n  p 

For  an  MOS-device  Laplace's  equation 

(2.7)  div(eVf)  -  0,  x  e  4 

holds  where  4  represents  the  oxide.  The  potential  and  the  electrical  displacement  are 
assumed  to  be  continuous  at  the  interface  3(lQg  -  3Q  O  34  and  electrons  and  holes  are  not 
allowed  to  penetrate  the  oxide  (n  and  p  vanish  in  4).  He  assume  that  the  device  has 
r  contacts  C|,...,Cr.  At  these  contacts  the  boundary  conditions 

(2.8)  4lc  -  ♦  +  UA 

are  prescribed  where  4  denotes  the  built  in  potential  at  the  contact  if  is  Ohmic, 

the  negative  flat-band  voltage  if  CA  is  a  gate  contact  on  the  metal-semiconductor 
interface.  Hark  function  if  is  a  Schottky  contact.  is  the  potential  applied 

to  C^. 

Also  n  and  p  are  prescribed  at  Ohmic  and  Schottky  contacts. 

•Bie  remaining  part  of  the  boundary  of  the  device  is  assumed  to  be  insulating,  that 
means  the  derivatives  of  4,  n  and  p  in  normal  direction  to  the  boundary  vanish  there. 


-V.VA  7T  V.1.'  u.’ y  r-r-T^r 


The  outflow  current  density  of  Cj  is  given  by 

(2.9)  I#  -  )  J  •  nd s  , 

H 

where  n  denotes  the  exterior  unit  normal  vector  to  . 

For  the  following  we  write  the  problem  (2.1)— (2.3)  (end  (2.7)  for  a  MOS  device) 
subject  to  the  mentioned  boundary  conditions  (and  interface  condition  for  a  MOS  device)  or 
a  scaled  version  of  it  in  operator  form: 

(2.10)  A(XfU^f...rO^)  ■  0|  z  m  (f ,n,p) 
while 

(2.11)  A  :  Z  *  *r  ♦  Y 

where  Z  and  Y  are  appropriate  (Banach)  spaces  of  functions.  In  the  sequel  we  are 
interested  in  the  determination  of  the  outflow  current  for  varying  contact  voltage 

0j  e  R  and  fixed  voltages  .  ,Uj_ ^,0^  j, . . . ,0r .  For  notational  simplicity  we  set 

0  i"  0^,  I  i“  ii#  we  assume  that  the  resulting  current-voltage  (I-u) -characteristic  is  a 
smooth  curve  in  R?  which  can  be  parameterized  as  follows 

(2.12)  O  -  0(s) ,  I  -  I(s) 
where 

(2.13)  i^(s)  +  i2(s)  -  1  (•  denotes  |^)  . 

Then  s  -  s0  is  the  length  of  the  arc  connecting  (0(sg),I(s0) )  and  (0(s),I(b)). 

The  goal  of  the  analysis  following  this  section  is  to  facilitate  the  numerical 
computation  of  the  curve  (2.12).  For  further  notational  simplicity  we  denote 

(2.14)  B(s,0)  -  A(s,01f ...,0j,1,0,Oj+1,...,0r> 

(all  contact  voltages  except  Oj  ■  0  are  fixed)  where 

(2.15)  B  i  Z  x  r  Y  . 
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3.  Parameterization 

The  arc length  parameterization  (2.13)  is  not  suited  for  computational  purposes  and 
therefore  we  use  the  approximation 

(3.1)  i(T0)(I  ‘  I(T0)>  +  “(To)tu  ‘  u(T0n  “  <T  "  V  *  0 

for  t  sufficiently  close  to  tq  assuming  that  (Uo>10)  !“  (u<Tq).I(Tq))  and  the  unit 

tangent-vector  (Uq,Iq)  i“  (U(t0), t(TQ) )  are  known. 

Noting  that  I  as  given  by  (2.9)  (for  t  «  i)  is  a  functional  of  z,  i.e. 


X  ■  l[z]  we  denote  (3.1)  by 


N(z,U,T)  -0,  N  s  Z  *  R  ♦  R 


and  augment  B(z,D)  *  0  by  (3.2),  that  means  we  solve 


P(y,t)  - 


N(y,T). 


0,  y  »  (z,U) 


where  T  is  in  a  sufficiently  small  neighborhood  of  tq.  He  will  show  that  the  I-U- 
characteristic  can  be  computed  by  solving  (3.3)  by  continuation  in  t  (at  least  locally 
about  t0)  even  if  0„  is  a  snap-back  voltage. 

A  geometrical  algorithm  for  the  determination  of  (U(t),I(t))  as  defined  by  (3.1)  is 
as  follows.  At  first  the  point  S0  on  the  tangent  to  (Oq,Xq),  whose  distance  to 
(O0,Ig)  is  It  -  Tg |  (if  T  -  Tq  >  0  the  vector  pointing  from  (U0,I0)  to  S0  is 
oriented  as  the  tangent  vector  (Uq,Xq)  and  if  t  -  tq  <  0  it  is  oriented  in  the 
opposite  direction)  is  determined.  The  normal  to  the  tangent  through  S0  is  intersected 
with  the  characteristic  curve  giving  <U(t),I(t)).  This  is  illustrated  in  Figure  2  at  a 
'regular  point*  (a  non-snap-back-point)  of  the  characteristic. 

Figure  3  demonstrates  the  reparameterization  at  a  snap-back-point  (U0,I0).  It  is 
intuitively  clear  that  (U(t),I(t))  can  be  determined  in  a  locally  unique  way  if 
It  —  Tq I  is  sufficiently  small. 

To  prove  this  we  define  'regular  points'  and  'simple  limit  points'  (i.e.  snap-back 
points)  mathematically  (see  H.  Keller  (1977)). 
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For  th«  reader  who  is  not  familiar  with  the  concepts  of  nonlinear  analysis  in  Banach 
spaces  we  remark  that,  for  understanding  the  implications  of  the  following,  it  suffices  to 
regard  B  as  a  function  from  RIC+1  into  R*  and  N  as  a  function  from  RK+2  into  R 
for  some  integer  K.  Then  Frechet  derivatives  (linearizations)  are  then  Jacobi  matrices 
of  corresponding  dimensions-  This  corresponds  to  analyzing  an  appropriate  discretization 
of  the  semiconductor  device  equations  (subject  to  the  mentioned  boundary  conditions)  and 
of  (3.1). 

We  define: 

1)  The  point  (z0,Ug)(zg  :*  z(sQ)  is  the  solution  of  B(z,U0)  -  0  with  Ig  -  I(z0))  is 

dB 

called  regular  solution  point  of  B(z,U)  -  0  if  the  Frechet  derivative  —  (z0,U0)  is 
one-to-one  and  onto  (i.e.  invertible). 

2)  The  point  (z0,Ug)  is  called  normal  limit  solution  point  of  B(z,U)  «0  if 
N(ff  (zq,Uq ) )  is  one  dimensional  (H( A)  denotes  the  null  space  of  the  operator  A), 

the  codimension  of  R(|^  (z0,Ug))  is  one  (R(A)  denotes  the  range  of  the  operator  A) 
and  if  ||  (z„,U0)  t  R(ff  (zg.Ug)). 

The  following  Theorem  is  along  the  lines  of  Keller's  (1977)  Theorem  3.3. 

Theorem  3.1:  Let  (z0,U0)  be  either  a  regular  solution  or  a  normal  limit  solution  of 

B(z,U)  -  0.  If  B  is  sufficiently  smooth  then  for  (T  -  tq)  sufficiently  small  there 

exists  a  unique  smooth  arc  of  solutions  (z( t),0(t)>  of  B(z,U)  -  0  passing  through 

3P 

(Z0,U0)  and  fulfilling  (3.3).  The  Frechet  derivatives  u}  are  one-to-one  and  onto 

along  this  arc. 


C0  3(z,u)  ( WV 


3 z  (z0'u0>  3u  <*0«V 


I(T0>  Tz  0,T0> 


o  o 

holds.  If  ( Zq , Bg )  is  a  regular  solution  point  of  B(z,U)  “  0,  then  —  (Zg,U0)  is 
one-to-one  and  onto.  To  show  that  (*g«ug»Tg)  is  one-to-one  and  onto  it  suffices 


to  prove  that 


The  invertibility  of  CQ  and  the  implicit  function  theorem  (see  Schwartz  (1969)) 


imply  Theorem  3.1.  □ 

The  Theorem  implies  that  -  at  least  for  It  -  t q I  sufficiently  small  -  the  solution 
arc  of  B(z,U)  m  0  can  be  computed  by  solving  (3.3)  (which  is  the  'reparameterized 
version  of  B(z,U)  »  0)  for  z  and  U  (by  continuation  in  t)  using  Newton's  method 
which  is  quadratically  convergent  along  the  whole  solution  arc. 


.  •  -  •  *.•  •  •  -  * 
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4.  Implementation 

For  the  numerical  eolution  of  the  fundamental  semiconductor  device  equations 
B(z,U)  -  0  the  problem  has  to  be  discretized  appropriately,  that  means  Bh(sh,U)  -  0  has 
to  be  solved  instead  of  B(z,U)  -  0  (h  represents  the  grid  parameter).  For  a  collection 
of  results  on  discretization  methods  see  Mock  (1983). 

He  now  assume  that  B^  has  already  been  chosen  and  that  Bh  :  *lc+1  ♦  H*.  Also  an 
appropriate  numerical  integration  rule  I^tz^]  has  to  be  used  to  discretize 

I[z]  «  J  I  •  nds. 

C. 

To  get  the  continuation  started  at  S  •  Sg  we  solve  B^(z^,Ug)  ”  0  obtaining  Zj^ 
(we  assume  that  (z^.Ug)  -s  (zh(T0),O(TQ))  is  a  regular  solution  of  B(zh,U)  -  0)  and 
(Ih  t«  lh [zh  )•  An  approximation  of  the  tangent  vector  to  the  I  -  0  characteristic  at 
(U0,I0)  can  be  obtained  by  solving  Bh**h'W0*  -  0  (with  AU0  “  U0  8utf*  8Ball) 

for  *h  ”  and  by  setting 


U0  "  uo 


with  <>»  -  /(UQ  -  Ug)2  +  (Iy, f*h01  "  rhf*hon 
Assume  now  that  we  already  solved 


V\5  "  IHt*h01 


VV'VV’  =“  'dr  VV'  d?  n '  K 1 


V*vV 

fh,T)  5 

N(zh,Oh,T) 

W 

are  known). 

1  from 

(4.2)  dT  Nh(zh<Uh,TK) 


5T  ‘VV'VV  *r  ‘WW  [Vv' 


*‘V  5^  ‘W  6,V 
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(4.4) (b) 


(QMVl).iMVln 
V(Vl))2  +  d*(Vl)2) 


h  way  to  determine  whan  an  update  is  desirable  is  to  require  that  the  length  of  the 
tangent  vector  (U* (Tr) ,1' (tr) )  differs  from  1  by  less  than  a  certain  prescribed  error 
tolerance. 


Of  course,  the  sign  of  the  increments  Atk  has  to  be  constant  throughout  the  whole 
continuation  process  (and  equal  to  the  sign  of  &UQ ) .  For  the  thyristor  characteristic 
with  Ug  »  Iq  «  0,  Tq  »  0  positive  increments  (AUg  >  0,Atk  >  0)  imply  that  one  'walks 
up'  the  forward  characteristic  (0  >  0)  and  negative  increments  (A(Jg  <  0,Atr  <  0) 
result  In  the  computation  of  the  reverse  characteristic  (U  <  0). 


The  increments  iTR  have  to  be  chosen  such  that  the  initial  guess  (4.3)  lies  in  the 
domain  of  attraction  of  the  iteration  method  for  (4.1)  at  t  -  tr+1.  For  a  sufficiently 
smooth  solution  arc 


(4.5) 


VW 

■T 

<tjc+i * 

<*V* 

w 

°h(Vl> 

C 

'Vi1 

2 

w 

'  *  lTK'TK+1 1 


holds.  Therefore  one  has  to  require  that 


(4.6) 


,4v 4 


where  represents  the  convergence  radius  of  the  iterative  method  for  the  approximate 

solution  of  (4.1)  at  t  m  t_+1.  Computable  estimates  for  4R+1  and  strategies  slightly 
different  froai  (4.6)  can  be  found  in  Den  Heiser  and  Rheinboldt  (1981). 

Assusm  now  that  a  Newton-type  procedure  is  used  as  iteration-method.  Then  a  linear 
system  of  the  form 


(4.7) 


has  to  be  solved  at  every  iteration  step  (and  also  for  solving  (4.2))  with 


Ah  *h 

«zh 

>* 

ch  d . 

.  «uh. 

MhJ 

(4.8) 


ffh 

^  “  3*h 
*Bh 


azh 

3U 


evaluated  at  some  iterate 

,  (t>  ..(i),  .  .  _  , 

(*h  i Uh  )  and  T  “  t 
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If  a  Gauss-solver  for  sparse  matrices  is  used  to  solve  systems  of  the  form  A^z^  =  rh 

(which  occur  In  solving  the  fundamental  semiconductor  device  equations  (2.10)  for 

prescribed  '■^-t-voltages)  then  it  can  easily  be  modified  to  solve  systems  of  the  form 

bh 

(5.7).  Tbe  K  +  1-st  row  )  and  the  K  +  1-st  column  (c^.d)  do  not  create 
additional  fill-in  during  the  elimination  process. 

If  a  band-solver  is  used  to  solve  systems  with  coefficient  matrices  Ah  then  systems 
of  the  form  (4.7)  can  be  solved  as  suggested  in  Keller  (1977),  that  is  by  solving  two 
systems  with  coefficient  matrix  simultaneously.  This,  of  course,  only  gives  accurate 

results  when  has  a  moderate  condition  number,  i.e.  when  Tj,  does  not  correspond  to  a 

limit  point.  However,  it  is  possible  to  continue  beyond  limit  points  when  using  the 
parameterization  (3.1)  (see  Keller  (1977),  Theorem  4.4). 


5.  A  Teat  Problem 


As  a  test  problem  we  applied  the  arclangth-continuation  procedure  to  the  one- 
dimensional  semiconductor  device  equations  modelling  a  silicon  pn-junctlon  in  the 

avalanche  case. 

For  simplicity  we  took  an  odd,  piecewise  continuous  doping  profile  (see  Figure  4). 


The  length  of  the  pn-junctlon  is  2t  m  5  x  10"3cm. 

The  avalanche  phenomenon  was  modelled  by  the  generation  rate 

(5.1)  R  -  “*<♦*)(  JJnl  +  Upl) 

with  Jn»Jp  given  by  (2.4),  (2.5)  resp.  The  following  expression  was  used  for  the  hole 
and  electron  ion  isolation  rate: 

o2 

(5.2)  a(*'(x))  -  (^expt-  |’j'.(x)'j')  >  <V®2  *  0  * 

The  electron  and  hole  mobilities  were  taken  equal  (also  the  electron  and  hole 
diffusivities).  The  obtained  results  should  not  be  regarded  as  physically  significant, 
they  however  illustrate  very  well  the  power  of  the  reparameteriration  technique. 
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We  remark  that  in  the  one-dimensional  case  the  current  is  constant  and  given  by 


(5.3)  1  5  q(0n(n'  -  p* )  -  Un(n  ♦  p)t*>  (for  Dn  »  Dp,  pn  -  up)  • 

A  scaled  version  of  (2.1)-(2.3)  (see  Markovich  (1983a)  for  the  scaling)  was  discretized  by 
using  the  three-point-scheme  for  Poisson's  equation  and  the  Scharfetter-Gummel  scheme  for 
the  continuity  equations  (see  Markovich  et  al.  (1983)). 

The  functional  Ih(zhl  was  obtained  by  discretizing  the  current  (5.3)  at  the  largest 
grid  point  xA  with  x^  <  1  also  using  the  Scharfetter-Gummel  discretization. 

Figure  5  shows  the  obtained  I-U-characteristic  for  -5UT  <  u  <  0  (0T  is  the  thermal 
voltage).  The  onset  of  avalanche  generation  occurs  at  U  ■  -1.5UT  and  the  current 
increases  linearly  (in  absolute  value)  beyond  that.  This  linear  increase  continues  as  far 
as  U  ■  -200Ut.  Then  the  increase  gets  faster  and  the  snap-back  occurs  at  U  «  -240OT 
(see  Figure  8).  The  continuation  of  the  solution  arc  beyond  the  limit  point  was  no 
problem. 

1 

Figures  7  and  8  show  J  a(+'(a))ds  over  the  applied  bias.  For  U  “  0  (equilibrium 
1  -1 

solution)  J  a(9')ds  <  1,  then  it  increases  until  it  equals  1.  This  happens  around  that 
-1 

O-value  at  which  avalanche  generation  starts,  i.e.  U  *  -1.5(1^  (compare  to  Pigure  5). 

Then  J  a(t')ds  remains  almost  constant  until  close  to  the  snap-back  voltage  (see  Figure 
-1 

4).  Slightly  before  the  snap-back  voltage  it  increases  and  gets  significantly  larger  than 

1.  It  was  often  claimed  (see,  for  example,  Sze  (1981))  that  'breakdown*  happens  at  that 
1 

voltage  at  which  J  a(9‘ Ids  reaches  one.  This  is  not  true  in  the  mathematical  sense  for 
-1 

our  staple  model,  there  is  a  continuous  branch  of  solutions  (of  the  one-dimensional 
semiconductor  problem  in  the  case  of  avalanche  generation)  which  emerges  at  the 
equilibrium  solution  (for  0  »  0)  and  which  contains  at  least  one  solution  for  every 
0  <  0  (see  Markovich  (1983b)).  The  situation,  however,  might  change  if  temperature  is 
introduced  as  unknown  quantity. 
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